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Abstract 



We propose a multi-scale stochastic volatility model in which a fast 
mean-reverting factor of volatility is built on top of the Heston stochastic 
volatility model. A singular pertubative expansion is then used to ob- 
Ch ' tain an approximation for European option prices. The resulting pricing 

C^_^ , formulas are semi-analytic, in the sense that they can be expressed as 

integrals. DifBculties associated with the numerical evaluation of these 

integrals are discussed, and techniques for avoiding these difficulties are 

provided. Overall, it is shown that computational complexity for our 

04 ■ model is comparable to the case of a pure Heston model, but our correc- 

^ , tion brings significant flexibility in terms of fitting to the implied volatility 

>D . surface. This is illustrated numerically and with option data. 

cn_ 

'^. ■ 1 Introduction 

^P . Since its publication in 1993, the Heston model [T^] has received considerable 

attention from academics and practitioners alike. The Heston model belongs 
to a class of models known as stochastic volatility models. Such models relax 
the assumption of constant volatility in the stock price process, and instead, 
allow volatility to evolve stochastically through time. As a result, stochastic 
^ ■ volatility models are able to capture some of the well-known features of the 

5^ , implied volatility surface, such as the volatility smile and skew (slope at the 

money). Among stochastic volatility models, the Heston model enjoys wide 
popularity because it provides an explicit, easy-to-computc, integral formula for 
calculating European option prices. In terms of the computational resources 
needed to calibrate a model to market data, the existence of such a formula 
makes the Heston model extremely efficient compared to models that rely on 
Monte Carlo techniques for computation and calibration. 
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Yet, despite its success, the Heston model has a number of documented short- 
comings. For example, it has been statistically verified that the model misprices 
far in-the-money and out-of-the- money European options [5] , |21j . In addition, 
the model is unable to simultaneously fit implied volatility levels across the 
full spectrum of option expirations available on the market [TU]. In particular, 
the Heston model has difficulty fitting implied volatility levels for options with 
short expirations [11]. In fact, such problems are not limited to the Heston 
model. Any stochastic volatility model in which the volatility is modeled as 
a one-factor diffusion (as is the case in the Heston model) has trouble fitting 
implied volatility levels across all strikes and maturities |llj . 

One possible explanation for why such models are unable to fit the implied 
volatility surface is that a single factor of volatility, running on a single time 
scale, is simply not sufficient for describing the dynamics of the volatility process. 
Indeed, the existence of several stochastic volatility factors running on different 
time scales has been well-documented in literature that uses empirical return 
data [J, [2], [3], [S], [5], [S], US], [13, [IH]- Such evidence has led to the 
development of multi-scale stochastic volatility models, in which instantaneous 
volatility levels are controlled by multiple diffusions running of different time 
scales (see, for example, [7]). We see value in this line of reasoning and thus, 
develop our model accordingly. 

Multi-scale stochastic volatility models represent a struggle between two 
opposing forces. On one hand, adding a second factor of volatility can greatly 
improve a model's fit to the implied volatility surface of the market. On the other 
hand, adding a second factor of volatility often results in the loss of some, if not 
all, analytic tractability. Thus, in developing a multi-scale stochastic volatility 
model, one seeks to model market dynamics as accurately as possible, while at 
the same time retaining a certain level of analyticity. Because the Heston model 
provides explicit integral formulas for calculating European option prices, it is 
an ideal template on which to build a multi-scale model and accomplish this 
delicate balancing act. 

In this paper, we show one way to bring the Heston model into the realm of 
multi-scale stochastic volatility models without sacrificing analytic tractability. 
Specifically, we add a fast mean-reverting component of volatility on top of the 
Cox-IngersoU-Ross (CIR) process that drives the volatility in the Heston model. 
Using the multi-scale model, we perform a singular perturbation expansion, as 
outlined in [7] , in order to obtain a correction to the Heston price of a European 
option. This correction is easy to implement, as it has an integral representation 
that is quite similar to that of the European option pricing formula produced 
by the Heston model. 

The paper is organized as follows. In Section [2] we introduce the multi-scale 
stochastic volatility model and we derive the resulting pricing partial differential 
equation (PDE) and boundary condition for the European option pricing prob- 
lem. In Section [3] we use a singular perturbative expansion to derive a PDE for 
a correction to the Heston price of a European option and in Section|4]we obtain 
a solution for this PDE. A proof of the accuracy of the pricing approximation is 
provided in Section [5] In Section [6] we examine how the implied volatility sur- 



face, as obtained from the multi-scale model, compares with that of the Heston 
model, and in Section [7] we present an example of calibration to market data. In 
Appendix |X] wc review the dynamics of the Heston Stochastic volatility model 
under the risk-neutral measure, and present the pricing formula for European 
options. An explicit formula for the correction is given in Appendix [B] and 
the issues associated with numerically evaluating the integrals-representations 
of option prices obtained from the multi-scale model are explored in Appendix 

El 

2 Multi-Scale Model and Pricing PDE 

Consider the price Xt of an asset (stock, index, ...) whose dynamics under the 
pricing risk-neutral measure is described by the following system of stochastic 
differential equations: 

dXt = rXtdt + l^tXtdWt, (2.1) 

St = ^tfiYt), (2.2) 

dYt = —{m-Yt)dt + v^/2\—dWf, (2.3) 

(- V e 



dZt = ti{e- Zt)dt + a^/ZtdW^. (2.4) 

Here, Wjf , VFf and W^ are one-dimensional Brownian motions with the corre- 
lation structure 

d{W\Wy)^ = p^ydt, (2.5) 

d{W^,W')t = p^.dt, (2.6) 

d{Wy,W'), = py^dt, (2.7) 

where the correlation coefficients pxy Pxz and py^ arc constants satisfying p^, < 
l.pi^ < I.pIz < 1, and ply + pl^ + p^^ - 2pxyPxzPyz < 1 in order to ensure 
positive definiteness of the covariance matrix of the three Brownian motions. 

As it should be, in (|2.ip the stock price discounted by the risk-free rate r is 
a martingale under the pricing risk neutral measure. The volatility Et is driven 
by two processes Yt and Zt, through the product y/Z^ f{Yt). The process Zt is a 
Cox-IngersoU-Ross (CIR) process with long-run mean 9, rate of mean reversion 
K, and "CIR-volatility" a. We assume that k, 9 and a are positive, and that 
2k9 > (7^, which ensures that Zt > ai all times, under the condition Zq > 0. 

Note that given Zt , the process Yt in (|2.3p appears as an Ornstcin-Uhlenbcck 
(OU) process evolving on the time scale e/Zt, and with the invariant (or long- 
run) distribution J\f{m, ^/^). This way of "modulating" the rate of mean rever- 
sion of the process Yt by Zt has also been used in [4] in the context of interest 
rate modeling. 

Multiple time scales are incorporated in this model through the parameter 
e > 0, which is intended to be small, so that Yt is fast-reverting. 



We do not specify the precise form of f{y) which wih not play an essential role 
in the asymptotic results derived in this paper. However, in order to ensure St 
has the same behavior at zero and infinity as in the case of a pure Heston model, 
we assume there exist constants ci and C2 such that < ci < f(y) < C2 < oo 
for all y S R. Likewise, the particular choice of an OU-like process for Yt is 
not crucial in the analysis. The mean-reversion aspect (or ergodicity) is the 
important property. In fact, we could have chosen Yt to be a CIR-like process 
instead of an OU-like process without changing the nature of the correction to 
the Heston model presented in the paper. 

Here, we consider the unique strong solution to (|2.1H2.4p for a fixed param- 
eter e > 0. Existence and uniqueness is easily obtained by (i) using the classical 
existence and uniqueness result for the CIR process Zt defined by (|2.4[) . (ii) 
using the representation (|5.18[) of the process Yt to derive moments for a fixed 
e > 0, (iii) using the exponential formula for Xt'. 
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We note that if one chooses f{y) = 1, the multi-scale model becomes e- 
independent and reduces to the pure Heston model expressed under the risk- 
neutral measure with stock price Xt and stochastic variance Zt: 

dXt = rXtdt + ^tXtdW^ , 
dZt = K{e- Zt)dt + a^/z'tdWt. 
d{W\W^)t = p^.dt. 

Thus, the multi-scale model can be thought of as a Heston- like model with a 
fast-varying factor of volatility, /(Yt), build on top of the CIR process Zt, which 
drives the volatility in the Heston Model. 

We consider a European option expiring at time T > t with payoff h{XT)- 
As the dynamics of the stock in the multi-scale model are specified under the 
risk- neutral measure, the price of the option, denoted by Ft, can be expressed 
as an expectation of the option payoff, discounted at the risk-free rate: 



Ft =E 



e-^-'-^-'^hiXr) Xt,Yt,Zt ^■. P'{t,Xt,Yt, Z, 



where we have used the Markov property of (Xt, Yt,Zt), and defined the pric- 
ing function P'^{t,x,y, z), the superscript e denoting the dependence on the 
small parameter e. Using the Feynman-Kac formula, P'^{t,x,y,z) satisfies the 
following PDE and boundary condition: 

C'P'{t,x,y,z) = 0, (2.8) 

^" = Q-^+^{x.Y.z)-r, (2.9) 

P'{T,x,y,z) = Hx), (2.10) 



where the operator C(^x,y.z) is the mfinitesimal generator of the process (Xj, 1^, Zt): 



d 1 ,2, X 2 c*^ -^ ^ d^ 



>n ^9 1 2 92 

z f. .d ^d^ 

+ 7(^(^-2^)5^ + - a^ 

It will be convenient to separate £'^ into groups of like-powers of l/\/e. To this 
end, we define the operators £0, ^1 and £2 as follows: 

^0 - -^|^ + (™-^)|^ (2-11) 

52 52 

£1 := Py^CTi/\/2— — - +pj;y!^V2/(y)x— — -, (2.12) 

With these definitions, C^ is expressed as: 

C = -Co + ^Ci+C2. (2.14) 

Note that Cq is the infinitesimal generator of an OU process with unit rate 
of mean-reversion, and £2 is the pricing operator of the Heston model with 
volatility and correlation modulated by f{y). 

3 Asymptotic Analysis 

For a general function /, there is no analytic solution to the Cauchy problem 
(|2.8H2.10|) . Thus, we proceed with an asymptotic analysis as developed in [7]. 
Specifically, we perform a singular perturbation with respect to the small pa- 
rameter e, expanding our solution in powers of y/e 

P' = Pn + V^Pi+eP2 + ... . (3.1) 

We now plug (|3?T|) and p.l4p into ([2^ and ()2.10p . and collect terms of equal 
powers of ^/e. 



The Order 1/e Terms 

Collecting terms of order 1/e we have the following PDE: 

= zCoPo- (3.2) 

We sec from (|2.1ip that both terms in £o take derivatives with respeet to y. In 
fact, £o is an infinitesimal generator an consequently zero is an eigenvalue with 
constant eigenfunctions. Thus, we seek Po of the form 

Po = Pa{t,x,z), 

so that (13.21) is satisfied. 



The Order 1/ ^/e Terms 

Collecting terms of order l/-y/e leads to the following PDE 

= zCoPi + zCiPo 

= zC^Pi. (3.3) 

Note that we have used that £iPo = 0, since both terms in Ci take derivatives 
with respect to y and Pq is independent oi y. As above, we seek Pi of the form 

Pi = Pi{t,x,z), 

so that (13.31) is satisfied. 



The Order 1 Terms 

Matching terms of order 1 leads to the following PDE and boundary condition: 

= zCaP2 + zCiPi+C2Po 

= zCaP2+C2Po (3.4) 

h{x) = Po{T,x,z). (3.5) 

In deriving (|3.4p we have used that £iPi = 0, since Ci takes derivative with 
respect to y and Pi is independent of y. 

Note that p.4|) is a Poisson equation in y with respect to the infinitesimal 
generator Cq and with source term C2P0', in solving this equation, (t,x,z) are 
fixed parameters. In order for this equation to admit solutions with reasonable 
growth at infinity (polynomial growth), we impose that the source term satisfies 
the following centering condition: 

= {C2P0) = {C2) Po, (3.6) 



where we have used the notation 

{g) := j 9{y)Hy)dy, (3.7) 

here <I> denotes the density of the invariant distribution of the process Yt , which 
we remind the reader is N{m, v^). Note that in (|3.6|) . we have pulled Po{t, x, z) 
out of the linear (•) operator since it does not depend on y. 

Note that the PDE (|3.6p and the boundary condition (|3.5p jointly define a 
Cauchy problem that PQ{t,x,z) must satisfy. 

Using equation p.4[) and the centering condition p.6p we deduce: 

P2 = --C-^{C2-{C2))Po, (3.8) 

z 

where £q is the inverse operator of Cq acting on the centered functions. 

The Order i/e Terms 

Collecting terms of order y^, we obtain the following PDE and boundary con- 
dition: 

= zLQP:i + zCiP2+C2Pi, (3.9) 

= Pi(r,x,2). (3.10) 

We note that P^{t^x^y,z) solves the Poisson equation (|3.9|) in y with respect 
to £9. Thus, we impose the corresponding centering condition on the source 
ZC1P2 + C2P1, leading to 

(£2)^1 = -{ZC1P2). (3.11) 

Plugging P2, given by (|3.8p . into equation p. lip gives: 

(£2)^1 = APo, (3.12) 

A := lzCi-^C^^{C2-{C2))Y (3.13) 

Note that the PDE (|3.12p and the zero boundary condition p.lOp define a 
Cauchy problem that Pi{t,x,z) must satisfy. 

Summary of the Key Results 

We summarize the key results of our asymptotic analysis. We have written the 
expansion p.ip for the solution of the PDE problem (|2.8H2.10"| . Along the way, 
he have chosen solutions for Pg and Pi which are of the form Pq = Po(t,a;,z) 
and Pi = Pi(t,x, z). These choices lead us to conclude that Po(i,a;,z) and 



Pi{t,x,z) must satisfy the following Cauchy problems 



and 



{C2)Po - 0, (3.14) 

Po{T,x,z) = h{x), (3.15) 



{C2)Pi{t,x,z) = APo{t,x,z), (3.16) 

PiiT,x,z) = 0, (3.17) 



where 



and ^ is given by p.l3p . Recall that the bracket notation is defined in p.7p . 

4 Formulas for PQ{t,x,z) and Pi{t,x,z) 

In this section we use the results of our asymptotic calculations to find explicit 
solutions for Po{t, x, z) and Pi (i, x, z). 

4.1 Formula for PQ{t,x,z) 

Recall that Po(^i2:, z) satisfies a Cauchy problem defined by equations (j3.14p 
and dSHD. 

Without loss of generality, we normalize / so that (^p) = 1. Thus, we 
rewrite (£2) given by (|3.18p as follows: 

P ■■= P.zif)- (4.2) 

We note that p^ < I since (/) < ^/^) = 1. So, p can be thought of as an effec- 
tive correlation between the Brownian motions in the Heston model obtained in 
the limit e — > 0, where (£2) = ^Hi the pricing operator for European options 
as calculated in the Heston model. Thus, we see that PQ{t,x,z) =: PH{t,x,z) 
is the classical solution for the price of a European option as calculated in the 
Heston model with effective correlation p — p^z if) ■ 

The derivation of pricing formulas for the Heston model is given in Appendix 



[XI Here, we simply state the main result: 

PH[t,x,z) = e-"-^^ I e-''"'G{T,k,z)h{k)dk, (4.3) 

T(t) = T-t, (4.4) 

q{t,x) = r(T-i)+loga;, (4.5) 

h{k) = I e'^'ih{e'i)dq, (4.6) 



G(r,fc,z) = e'^Kfc)+^-D(r,fc)^ ^47) 

C(r,fc) = ^i^[^ + p^ka + d{k))T-2\og{^—^^^^P-y^,{A.%) 

D(rk) = !i±P!fc^±M f ^^£1!!!^^ (4 9) 

d(A:) = y/a^{k^-ik) + {K + pika)^, (4.10) 

We note that, for certain choices of h, the integral in (|4.6p may not converge. 
For example, a European call with strike K has /i(e'') = (e"* — /\ )+. In this case, 
the integral in (|4.6p converges only if we set k = fc,. + ifc^ where ki > 1. Hence, 
when evaluating (|4.31 14. 6|) one must impose k = kr + iki, kr > I and dk = dkr- 



4.2 Formula for Pi(t,x,z) 

Recall that Pi {t, x, z) satisfies a Cauchy problem defined by equations (j3.16p 
and (|3.17p . In order to find a solution for Pi(t, x, z) we must first identify the 
operator A. To this end, we introduce two functions, 4>{y) and '!/'(y)j which solve 
the following Poisson equations in y with respect to the operator Cq: 

A0 = \{f-{f)). (4.12) 

^oV- = /-(/)• (4.13) 

From equation p.l3p we have: 

A = IzCi-^C^^ {C2 - {C2)) 



IzLy-Lq^Px70z[j- {f))x 



32 



dxdz 



- ( '^^^^y)^"§^) + p^^'^' {^'^'^y^'^^) 



Using the definition (|2.12p of £i, one deduces the following expression for A: 

Q3 g3 



A = VlZx'^^-^^TT^ +V2ZX- 



dzdx^ dz^dx 

+ V3ZX— U^— - + ^42 — 

OX \ ax^ ) az 

V2 = p^zPyz<y'^V\/2{%j/) , 

1/3 = p^yv^/2{fcp'), 

Vi = pxyPxzCrv\f^{fi'') ■ 

Note that we have introduced four group parameters, Vi, i 

constants, and can be obtained by calibrating our model to the market as will 

be done in Section [T] 

Now that we have expressions for A, Ph, and Ch, we are in a position to 
solve for Pi{t,x,z), which is the solution to the Cauchy problem defined by 
equations p.l6p and p.l7p . We leave the details of the calculation to Appendix 



x^y 

dx ) 


(4.14) 




(4.15) 




(4.16) 




(4.17) 




(4.18) 


K,i = l. 


. . 4, which are 
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[Bl Here, we simply present the main result. 

Pi{t,x,z) = ^ f e-''"^(K0fo{T,k)+zMT,k)) 

xG{T,k,z)h{k)dk, (4.19) 

Tit) = T-t, 
q{t,x) = r(T -t)+ log x, 

foir,k) = / his,k)ds, (4.20) 



Jo 
7i(T,fc) = r h{s,k)e^^^'^''^ds, (4.21) 







C(r,fc) = ^^[{n + p^ko + d{k))r-2\og(^—^^^^ 

^, ,, K + p2fc(T + d(A:) / l-e^^'W 
D(T,fc) = 5 



1 -g{k)e^d{k)J ' 



./ , ^ . •, j.,NNl-9(fc), /g(fc)e"''(*'-)-l 

A(t, fc, s) = {k + paik + d{k)) ,.,x ./^ log 



d(%(fc) "^ VsWe"''^''') -1, 
+d(fc) (r - s) , (4.22) 

d{k) = ^/(T^{k^ - ik) + {k + pikaY , 

, , K + pika + d(k) 

^(^) ^ . + Ma-d(fc) ' 

6(T,fc) = -{VlD{T,k){-P+ik)+V2D^{T,k){-ik) 

+V3 {ik^ + k^) + ViD{T, k) (-fc2)) . (4.23) 

Once again, we note that, depending on the option payoff, evaluating equation 
(|4.19p may require setting k = kj. + iki and dk = dk^, as described at the end 
of subsection 14.11 

5 Accuracy of the Approximation 

In this section, we prove that the approximation P'^ ^ Pi^ + \f^P\-, where Pq 
and P\ are defined in the previous sections, is accurate to order e" for any given 
a £ (1/2, 1). Specifically, for a European option with a smooth bounded payoff, 
/i(x), and with bounded derivatives, we will show: 

|P^(t,x,y,z)-(Po(i,x,z) + ViPi(t,a;,z))| < Ce", (5.1) 

where C is a constant which depends on (y, z), but is independent of e. 
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We start by defining the remainder term R'^it, x, y, z): 

R" = {Po + V^Pi + e^2 + ^V^-Ps) - P'- (5.2) 

Recalling that 

= C'-P\ 

= zCnPn, 

= zCaPi + zCiPo, 

^ zCaP2 + zCiPi + C2P0, 

= zCaP:i + ZC1P2 + C2P1, 

and applying C^ to R"^, we obtain that i?^ must satisfy the following PDE: 

CR" = C- {Po + V^Pi + eP2 + eV^Pa) - C P' 

^Co + -j=Ci + C2J (Po + V^-Pi + eP2 + (VePi) 

= e {zCiPa + C2P2 + ,f€C2Pi) 

= eP^ (5.3) 

F' := zCiP^+C2P2 + V€C2Pi, (5.4) 

where we have defined the e-dependent source term F'^{t, x, y, z). Recalling that 

P''{T,x,y,z) = h{x), 
Po{T,x,z) ~ h{x), 
Pi{T,x,z) = 0, 

we deduce from (15.21) that 



R\T,x,y,z) = eP2{T,x,y,z) + ey/eP3{T,x,y,z) 

= eG'ix,y,z), (5.5) 

G''{x,y,z) := P2{T,x,y, z) + ^/eP3{T,x,y,z), (5.6) 

where we have defined the e-dependent boundary term G'^{x, y, z). 

Using the expression (|2.9p for £^ we find that R'^{t,x,y, z) satisfies the fol- 
lowing Cauchy problem with source: 

^^+Cx,y.z-Ar' = eF', (5.7) 

R'iT,x,y,z) - eG'ix,y,z). (5.8) 



f2 



Therefore R^ admits the foUowing probabilistic representation: 



R'(t,x,y,z) = eE 



e-^iT-t)G^XT,YT,ZT) 






(5.9) 



In order to bound R'^(T, x, y, z), we need bounds on the growth of F'^{t, x, y, z) 
and G'^{x,y,z). From equation (|5.6p we see that G'^{x,y,z) contains the func- 
tions P2{t, X, y, z) and P3(i, x, y, z). And from equation (|5.4I) we see that F'^{t, x, y, z) 
contains terms with the hnear operators, £i and £2, acting on P2{t,x,y, z) 
and F3(t,x,y, z). Thus, to bound F^(t,x,y, z) and G'^{x,y, z), we need to ob- 
tain growth estimates for P2{t,x,y,z), Pz{t,x,y,z) and growth estimates for 
P2{t, X, y, z) and Psit, x, y, z) when hnear operators act upon them. To do this, 
we use the following classical result, which can be found in Chapter 5 of [3- 

Lemma 5.1. Suppose Cqx = 5; (.9) ~ and \g{y)\ < Ci(l + |y|"), then \x{y)\ < 
C2(l + I2/I") for some G2. When n = we have \x{y)\ < C2(l + log(l + |y|)). 

Now, by continuing the asymptotic analysis of Section|3l we find that P2 (t, x, y, z) 
and Psit, X, y, z) satisfy Poisson equations in y with respect to the operator, Cq. 
We have 

CoP2it, X, y,z) = - (-£2 + (£2)) Po{t, X, z), 

z 

CoPsit, X, y,z) = - {-C2 + (£2)) Pi {t, X, z) + {-CiP2(t, X, y, z) + {CiP2{t, .x, y, z))) . 

z 

Also note, for any operator, A4, of the form 

M = ^Y\x"^^^^-^^, (5.10) 

wehaveA^£o = £oA^ 7 because £0 does not contain x or z. Hence, A^P2(^, x, y, z) 
and MP3{t,x,y, z) satisfy the following Poisson equations in y with respect to 
the operator, Lq 

£0 {MP2{t, X, y, z))=M-^ (-£2 + (£2)) Po(i, x, z), (5.11) 

Co{MPi{t,x,y,z)) = M- {-C2 + {C2)) Pi{t,x,z) 

z 

+ M {-CiP2{t, X, y, z) + (£iP2(t, X, y, z))) . 
Let us bound functions of the form M.Po{t,x,z). Using equations (|4.3p and 
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(|5.10p , and recalling that q = rr + log x and G = q'^+'^d ^ ^g have 



MPq 



27r 



27r 



27r 



Af 






"(i) 






-ikq 



^C{TM,z)+zD(T,k,z) 



dz 



h{k)dk 



N n{j) 



e-*'=« II II {-ik -l + l)\ ({D{t, k, z)r e^("-'=-^)+^^("^'=^^)) h{k)dk 

N n(i) \ 

n n ("*^ - ^ + 1) (^(^' '^^ -2))" e-''^'«G(r, fc, z)h{k)dk. 



We note the following: 

• By assumption, the option payoff, h{e'^) G S, the Schwartz class of rapidly 
decreasing functions. It is a fact that the Fourier transform, h{k) € 5 as 



well. This implies that 



fc™/i(fc) 



< cxj for all integers, m. 



G{t, k, z) <l for all r e [0, T], fc e M, z e K+. This follows from the fact 
that G{t, fc, z) is the characteristic function, ¥\cyi^{ikQT)\Xt ~ x, Zt ~ z]. 

• There exists a constant, C, such that |£'(r, fc)| < C(l + |fc|) for all r g 

[o,r]. 



It follows that for any Ai of the form (|5.10p wc have the following bound on 
MPo{t,x,z) 



\MPo{t,x,z)\< 



2tt 



N n(j) 

niiHfc-^+i) 



< 



N n(j) 

nriHfc-^+i) 



\D{T,k)r\e-'^'i\ G{T,k,z) h{k) 



dk 



dk :— C < oo, 



(5.12) 



The constant C depends on A4, but is independent of {t,x,z). Using simi- 
lar techniques, a series of tedious but straightforward calculations leads to the 
following bounds 



\MPi{t,x,z)\ <C{l + z), 



MPoit,x,z) 
MPi(t,x,z) 



<Cil + z), 
<C(l + z2), 
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where, in each case, C is some finite constant which depends on M.^ but is 
independent of {t,x,z). We are now in a position to bound functions of the 
form A4P2{t, x, y, z) and A^P3(i, x, y, z). From equation (j5.1ip we have 

£0 {MP2{t, X, y, z))=M^ (-£2 + {C2)) Po{t, X, z) 

= \{-f{y) + {f))MiPa{t,x,z) 

+ p^.a i-fiy) + if)) M2Po{t,x, z) 
=■ 9{t,x,y,z), 



where Mi are of the form (|5.10p . Now using the fact that /(y) is bounded and 
using equation (|5.12p we have 

\9{t,x,y,z)\ < C, 

where C is a constant which is independent of {t,x,y,z). Hence, using lemma 
15.11 there exists a constant, C, such that 

\MP2{t,x,y,z)\<C{l+\og{l + \y\)). 

Similar, but more involved calculations, lead to the following bounds: 



\MP^{t,x,y,z 


1, 


^MP2it,x,y,z) 


<C(l + log(l + |y|))(l + z), 


(5.13) 




—MP2{t,x,y,z) 


<C, 


^MP2it,x,y,z) 


■} 


■^MP3it,x,y,z) 


<C(l + z), 




^^MP3{t,x,y,z) 


<C(l + log(l + |2/|))(l + z2), 




d 

dy 




-—MP3it,x,y,z) 
at 


<C(l + z2). (5.14) 



We can now bound G"^(a;, y, z). Using equation (|5.6p we have 

\G'{x,y,z)\ <\P2{T,x,y,z)\+ ^e\P3{T,x,y,z)\ 

< Ci(l + log(l + |y|)) + V^C2(1 + log(l + \y\)){l + z) 
<C(l + log(l + |y|))(l + z). (5.15) 

Likewise, using equation (j5.4[) . we have 

\F'{t,x,y,z)\<z\CiP3{t,x,y,z)\ + \C2P2{t,x,y,z)\+^e\C2P3{t,x,y,z)\. 
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Each of the above terms can be bounded using equations (|5. 13115.1^ . In partic- 
ular we find that there exists a constant, C, such that 

\F'{t,x,y,z)\ < C(l +log(l + |y|))(l + z^). (5.16) 

Using (|5.9p . the bounds (|5.15p and (|5.16p . Cauchy-Schwarz inequality, and mo- 
ments of the e- independent CIR process Zt (see for instance [U]), one obtains: 



\R'{t,x,y,z)\<eC{z) ll+Et,y,,\YT\+ J Et,j,,,|y,|ds , (5.17) 

where Kt,y,z denotes the expectation starting at time t from Yt ^ y and Zt — z 
under the dynamics (|2.3I) - (|2.4I) . Under this dynamics, starting at time zero 
from y, we have 

Yt^m + {y~ m)e-i fo z.ds ^ Y2^^-\ /„' z^du f ^\ f^ ^-^^v^.dW^,. 

(5.18) 

Using the bound established in Appendix [Cl we have that for any given a. € 
(1/2, 1) there is a constant C such that. 

E|rt|<Ce"-\ (5.19) 

and the error estimate (15.11) follows. 



Numerical Illustration for Call Options 

The result of accuracy above is established for smooth and bounded payoffs. 
The case of call options, important for implied volatilities and calibration de- 
scribed in the following sections, would require regularizing the payoff as was 
done in [^ in the Black-Scholes case with fast mean-reverting stochastic volatil- 
ity. Here, in the case of the multi-scale Heston model, we simply provide a 
numerical illustration of the accuracy of approximation. The full model price 
is computed by Monte Carlo simulation and the approximated price is given by 
the formula for the Heston price Pq given in Section 14.11 and our formulas for 
the correction ^/e P\ given in Section 14.21 Note that the group parameters Vi 
needed to compute the correction are calculated from the parameters of the full 
model. 

in Table [1] we summarize the results of a Monte Carlo simulation for a 
European call option. We use a standard Euler scheme, with a time step of 
10~^ years-which is short enough to ensure that Zt never becomes negative. 
For each value of e wc run 10^ sample paths. The parameters used in the 
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simulation are: 

X = 100, z = 0.24, r = 0.05, K=l,e = l,a = 0.39, p^^ = -0.35, 

y = 0.06, m = 0.06, v — 1, Pxy = —0.35, pyz = 0.35, 

T = l,K= 100, 

and /(y) = gV-"^-'^ go that (/^) = 1. Note that ahhough / is not bounded, it 
is a convenient ehoice because it allows for analytic calculation of the four group 
parameters Vi given by (j4.15l[4?T8l) . We only display the value of the largest 
one, Ve V3, which controls the correction of the skew due to the presence of pxy 
We note that the value of -y/e V3 calibrated to data from the S&P500 in Section 
[7|is even smaller than those displayed in the Table. 

The first line of Table [T] corresponds to the case of a pure Heston model 
(e = 0). Therefore, the value Pq = 21.0831 is exact (computed with analytic 
formulas), and it gives us a calibration of the empirical error due to the Monte 
Carlo simulation {umc = 0.1166). Note that this empirical error is consistent 
across the values of e used in the Table. 

As expected, the approximated price Pq + \f^P\i converges, as e — >■ 0, to the 
pure Heston price, and the approximation falls within one standard deviation 
of the Monte Carlo price for e < 10"'^. This illustrates the accuracy of our 
approximation for call options. 

e ^T^y-i Po + V^Pi Pmc ^mc \Po + yfePi-PMc\ 






0.0000 


21.0831 


21.1591 


0.1166 


0.0760 


lo-'' 


0.0096 


21.0055 


21.0045 


0.1153 


0.0010 


10-3 


0.0303 


20.8546 


20.7824 


0.1136 


0.0722 


10-2 


0.0959 


20.3752 


18.8250 


0.1015 


1.5502 


10-1 


0.3033 


18.8538 


14.8158 


0.0866 


4.0380 



Table 1: Results of a Monte Carlo simulation for a European call option. 



6 The Multi-Scale Implied Volatility Surface 

In this section, we explore how the implied volatility surface produced by our 
multi-scale model compares to that produced by the Heston model. To begin, 
we remind the reader that an approximation to the price of a European option 
in the multi-scale model can be obtained through the formula 



P" r 


- Po + V~^Pi 




= Ph+PI 


Pi ■■ 


= V~ePi, 
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where we have absorbed the %fk into the definition of P{ and used Pq = Pk , the 
Heston price. Form the formulas for the correction Pi, given in Section l42l it 
can be seen that Pi is hnear in V^, i = 1, • • • , 4. Therefore by setting 

V^ = yfeV, i = 1...4, 

the smaU correction PI is given by the same formulas as Pi with the Vi replaced 
by the Vf. 

It is important to note that, although adding a fast mean-reverting factor of 
volatility on top of the Heston model introduces five new parameters {v, m, e, 
Pxy, Pyz) plus an unknown function / to the dynamics of the stock (see (|2.2|) and 
(|2.3p ). neither knowledge of the values of these five parameters, nor the specific 
form of the function / is required to price options using our approximation. The 
effect of adding a fast mean-rcvcrting factor of volatility on top of the Heston 
model is entirely captured by the four group parameters V^ , which are constants 
that can be obtained by calibrating the multi-scale model to option prices on 
the market. 

By setting F/ = for z = 1, • • • , 4, we see that PI = 0, P^ = Ph, and the 
resulting implied volatility surface, obtained by inverting Black-Scholes formula, 
corresponds to the implied volatility surface produced by the Heston model. If 
we then vary a single V^ while holding V^ = Q for j ^ i, we can see exactly 
how the multi-scale implied volatility surface changes as a function of each of 
the VI . The results of this procedure are plotted in Figure [TJ 

Because they are on the order of -^/e, typical values of the V^ are quite small. 
However, in order to highlight their effect on the implied volatility surface, the 
range of values plotted for the V^ in Figure [T] was intentionally chosen to be 
large. It is clear from Figure [1] and from equation ()4.23|) that each V^ has a 
distinct effect on the implied volatility surface. Thus, the multi-scale model 
provides considerable flexibility when it comes to calibrating the model to the 
implied volatility surface produced by options on the market. 

7 Calibration 

Denote by 8 and <& the vectors of unobservable parameters in the Heston and 
Multicale approximation models respectively. 

e = {k,p,(7,9,z), 

$ = {K,p,a,9,z,Vi\V^\V,\Vi). 

Let a{Ti,Kj(j^)) be the implied volatility of a call option on the market with 
maturity date Tj and strike price ^j(i). Note that, for each maturity date, 
Tj, the set of available strikes, {Kj^^i)}, varies. Let (7H{Ti,Kj{^i),Q) be the 
implied volatility of a call option with maturity date Ti and strike price ii^j(i) as 
calculated in the Heston model using parameters 6. And let aM{Ti,Kj^i),<^) 
be the implied volatility of call option with maturity date Ti and strike price 
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(a) 



(b) 





(c) 



(d) 



Figure 1: Implied volatility curves are plotted as a function of the strike price 
for European calls in the multi-scale model. In this example the initial stock 
price is X = 100. The Heston parameters are set to z = 0.04, = 0.024, k = 3.4, 
a = 0.39, pxz = —0.64 and r = 0.0. In subfig ure|(a)| we vary only V^, fixing 
y/ = for i ^ 1. Likewise, in subfigures |(b)[ |(c)| and [(d)[ we vary only V2, 
V^ and V^ respectively, fixing all other Vf = 0. We remind the reader that, in 
all four plots, y/ = corresponds to the implied volatility curve of the Heston 
model. 



i^j(i) as calculated in the multi-scale approximation using parameters $. 

We formulate the calibration problem as a constrained, nonlinear, least- 
squares optimization. Define the objective functions as 

We consider Q* and $* to be optimal if they satisfy 

A2,(e*) = minAl,(e), 
Ai,(<I>*) = minAi,($). 

It is well-known that that the objective functions, A^ and Aj^j, may exhibit a 
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number of local minima. Therefore, if one uses a local gradient method to find 
0* and $* (as we do in this paper) , there is a danger of ending up in a local 
minima, rather than the global minimum. Therefore, it becomes important to 
make a good initial guess for Q and $, which can be done by visually tuning 
the Heston parameters to match the implied volatility surface and setting each 
of the Vf = 0. In this paper, we calibrate the Heston model first to find 6*. 
Then, for the multi-scale model we make an initial guess $ = (B*, 0, 0, 0, 0) (i.e. 
we set the Vf = and use 6* for the rest of the parameters of <J>). This is a 
logical calibration procedure because the F/, being of order -y/i, are intended to 
be small parameters. 

The data we consider consists of call options on the S&P500 index (SPX) 
taken from May 17, 2006. We limit our data set to options with maturities 
greater than 45 days, and with open interest greater than 100. We use the yield 
on the nominal 3-month, constant maturity, U.S. Government treasury bill as 
the risk-free interest rate. And we use a dividend yield on the S&P 500 index 
taken directly from the Standard & Poor's website (www.standardandpoors.com). 
In Figures [5] through [SJ wc plot the implied volatilities of call options on the 
market, as well as the calibrated implied volatility curves for the Heston and 
multi-scale models. Wc would like to emphasize that, although the plots are 
presented maturity by maturity, they are the result of a single calibration pro- 
cedure that uses the entire data set. 

From Figures [5] through |51 it is apparent to the naked eye that the multi- 
scale model represents a vast improvement over the Heston model-especially, 
for call options with the shortest maturities. In order to quantify this result we 
define marginal residual sum of squares 

where N{Ti) is the number of different calls in the data set that expire at time 
T, (i.e. N{Ti) = #{A'j(j)}). A comparison of Aj^iTi) and A^^(Ti) is given in 
Table [21 The table confirms what is apparent to the naked eye-namely, that 
the multi-scale model fits the market data significantly better than the Heston 
model for the two shortest maturities, as well as the longest maturity. 

A Heston Stochastic Volatility Model 

There are a number of excellent resources where one can read about the Hes- 
ton stochastic volatility model — so many, in fact, that a detailed review of the 
model would seem superfluous. However, in order to establish some notation, 
we will briefly review the dynamics of the Heston model here, as well as show 
our preferred method for solving the corresponding European option pricing 
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T^ 



- 1 (days) 

~65 
121 
212 
303 
394 
583 
947 



29.3 X 10"" 

10.2 X 10"^ 
4.06 X 10"^ 
3.93 X 10"6 
7.34 X 10-6 

11.3 X 10-6 
3.31 X 10-6 



ALm) 



Alm)/Al{T,) 



7.91 X lO"'* 
3.72 X 10-6 
8.11 X 10-6 
3.51 X 10-6 
5.17 X 10-6 
9.28 X 10-6 
1.47 X 10-6 



3.71 
2.73 
0.51 
1.12 
1.42 
1.22 
2.25 



Table 2: Residual sum of squares for the Heston and the Multi-Scale models at 
several maturities. 
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Figure 2: SPX Implied Volatilities from May 17, 2006 



problem. The notes from this section closely follow [20] . The reader should be 
aware that a number of the equations developed in this section are referred to 
throughout the main text of this paper. 

Let Xt be the price of a stock. And denote by r the risk-free rate of interest. 
Then, under the risk-neutral probability measure, P, the Heston model takes 
the following form: 



dXt 
dZt 



rXtdt + ^/Zt XtdW^ , 
K{e-Zt)dt + ay/z'tdWt 



pdt. 
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Figure 3: SPX Implied Volatilities from May 17, 2006 
Days to Maturity = 212 
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Figure 4: SPX Implied Volatilities from May 17, 2006 



Here, W'^ and W^ are one-dimensional Brownian motions with correlation p, 
such that IpI < 1. The process, Zt-, is the stochastic variance of the stock. And, 
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Figure 5: SPX Implied Volatilities from May 17, 2006 
Days to Maturity = 394 
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Figure 6: SPX Implied Volatilities from May 17, 2006 

K, 9 and a are positive constants satisfying 2k9 > cr^; assuming Zq > 0, this 
ensures that Zt remains positive for all t. 
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Figure 7: SPX Implied Volatilities from May 17, 2006 
Days to Maturity = 947 
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Figure 8: SPX Implied Volatilities from May 17, 2006 

Wc denote by Ph the price of a European option, as calculated under the 
Heston framework. As we arc already under the risk-neutral measure, we can 
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express Ph as an expectation of the option payoff, h{XT), discounted at the 
risk-free rate. 



PH{t,x,z) 



E 



'r(T-t) 



h{Xj 



Xf — X, Zt 



Using the Fcynman-Kac formula, wc find that PH{t,x, z) must satisfy the fol- 
lowing PDE and boundary condition: 



CnPHit^x^z) -- 


= 0, 


(A.l) 


PHiT,x,z) = 


= K^), 


(A.2) 


Ch -- 


d 9 1 2 52 

dt '+'^ax + 2^^ 9x2 
92 

A-nrrzx . 


(k.?,\ 



dxdz 

In order to find a solution for Pnit^x^z), it will be convenient to transform 
variables as follows: 



r(t) = T-t, 
q{t,x) = r(r-i)+logx, 
PH{t,x,z) = P^(T(t),g(i,.T),z)e- 



r(t) 



This transformation leads us to the following PDE and boundary condition for 
P'Hir,q,z): 



C'HP^{T,q,z) = 0, 



92 



d_ 1 / 92 Q 
dr 2 ydq"^ dqj dqdz 

1 2 5^ ^. ^ d 



+ 0^^ ^"a~2 + i^{G - z) 



d 



9z' 



(A.4) 



P'H{^,q,z) 



hie''). 



We will find a solution for P^ through the method of Green's functions. Denote 
by 5{q) the Dirac delta function, and let G{T,q, z), the Green's function, be the 
solution to the following Cauchy problem: 



C'HG{T,q,z) = 0, 
G(0,(7,z) = Siq). 



(A.5) 
(A.6) 



Then, 



Pkir,q,z) 



G{t, q — p, z)h{e^)dp. 
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Now, let Ph{t, k, z), G{t, k, z) and h{k) be the Fourier transforms of P'u{t, q, z) 
G{t, q, z) and /i(e') respectively. 

PH{T,k,z) = f e'''''P'H{T,q,z)dq, 

JR 

G{T,k,z) = f e''"iG{T,q,z)dq, 

JR 

h{k) = [ e''"^h{e'')dq. 

Jr 

Then, using the convolution property of Fourier transforms we have: 

P^{T,q,z) = ^ [ e-'''^PH{T,k,z)dk 

^'^ Jr 
= — f e-''"'G{T,k,z)h{k)dk. 

Multiplying equations (jA.5[) and ()A.6P by e**^* and integrating over M in g', we 
find that G{t, k, z) satisfies the following Cauchy problem: 

£ffG(T,fc,z) = 0, (A.7) 

+ [tiO - (k + paik) z) — , 

G(0,fc,z) = 1. (A.8) 

Now, an ansatz: suppose G{t, fc, z) can be written as follows: 

G{t, k, z) = eC^(^.fe)+^^(^^'^). (A.9) 

Substituting (jA.9[) into (jA.7[) and (jA.8[) . and collecting terms of like-powers of 
z, we find that C{t, k) and D{t, k) must satisfy the following ODE's 

dC 

-— (T,fc) = K9D{T,k), (A.IO) 

dr 
C(0,fc) = 0, (A.ll) 

^(T,fc) = ^a^D^{T,k)-{K + paik)D{T,k) + ^{-k'^ + ik),{A.l2) 
D{0,k) = 0. (A.13) 

Equations (|A.10|) . (JA.III) . (|A.12|) and (jA.13|) can be solved analytically. Their 
solutions, as well as the final solution to the European option pricing problem 
in the Heston framework, are given in (|4.3H4.11"|) . 
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B Detailed solution for Pi{t,x,z) 

In this section, we show how to solve for Pi(t, x, z), which is the solution to the 
Cauchy problem defined by equations p.l6p and p.l7p . For convenience, we 
repeat these equations here with the notation Ch ~ {C2) and Pr = Pa' 

£HPi{t,x,z) = APHit,x,z), (B.l) 

Pi{T,x,z) = 0. (B.2) 

We remind the reader that A is given by equation (|4.14l) , Ch is given by equation 
(|4.1|) . and Pf/(t, a;, z) is given by equation ()4.3p . It will be convenient in our 
analysis to make the following variable transformation: 

Piit,x,z) = P[{T{t),q{t,x),z)e-^\ (B.3) 

r{t) = T-t, 
q{t,x) ^ r{T - t)+ log X, 

We now substitute equations ()4.3p . ()4.14|) and (jB.SP into equations (jB.ip and 
(|B.2p . which leads us to the following PDE and boundary condition for P[{t, q, z): 

£'HPi{T,q,z) = A'^ f e-''"iG{T,k,z)h{k)dk, (B.4) 

27r J 

d 1 f d^ d\ ^^ 



dr 2 ydq"^ dqj dqdz 

1 a^ 

A' = Viz— { TT^ - ^] + V2Z- 



dz \dq'^ dqJ dz'^dq 

QZ Q2 X QZ 



\dq^ dq^ ) dzdq^' 

P[iO,q,z) - 0. (B.5) 

Now, let Pi{t, k, z) be the Fourier transform of P[{t, q, z) 

Pl{r,k,z) = f e'"^Pi{T,q,z)dq. 

Jr 

Then, 

Pi{r,q,z) = ^ [ e-''"^Pi{T,k,z)dk. (B.6) 
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Multiplying equations (jB.4p and (jB.Sp by e**^' and integrating in q' over R, we 
find that Pi(r, k, z) satisfies the following Cauchy problem: 

CHPi{T,k,z) = AG{T,k,z)h{k), (B.7) 



d 1 . .9 ..N 1 , a' 



2 , „-l\ , -^,2 



,2 



/:if = -j- + -z{-k^ + ik) + 



a z- 



dr 2 ^ '2 9z2 

d 

+ {k9 - (k + paik) z) — , 

\ oz az^ 

Pi(0,fc,z) = 0. (B.8) 

Now, an ansatz: we suppose that Pi(t, fc, z) can be written as 

Pi(r, fc, z) = (^efoiT, k) + zfi{T, fc))) GiT, k, z)h{k). (B.9) 

We substitute (jB.9p into (jB.7p and (jB.Sp . After a good deal of algebra (and in 
particular, making use of (jA.lOp and (|A.12p '). we find that /o(t, k) and /i(t, k) 
satisfy the following system of ODE's: 

^(T,fc) = a(T,fc)/i(T,fc)+6(r,fc), (B.IO) 

7i(0,fc) = 0, (B.ll) 

^(r,fc) = A(r,fc), (B.12) 

7o(0,fc) = 0, (B.13) 

a{T, k) = a D{t, k) — {k + paik) , 

&(T,fc) = -{VlD{T,k){-~k^ +ik)+V2D'^{T,k){-ik) 

+^3 (ifc3 + e) + F4Z?(t, fc) (-fc2)) , 



where I?(t, fc) is given by equation 

Equations (JB.10| - |BT13|) can be solved analytically (to the extent that their 
solutions can be written down in integral form). The solutions for /o(t, k) and 
/i(t, fc), along with the final solution for Pi(i, cc, z), are given by (|4.19H4.251) . 

C Moment Estimate for Yt 

In this section we will derive a moment estimate for Yj, whose dynamics under 
the pricing measure are given by equations (|2.3[ 12.41 12. 7p . Specifically, we will 
show that for all a £ (1/2, 1) there exists a constant, C (which depends on a 
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but independent of e), such that E|l^| < Ce"~^. 

We will begin by establishing some notation. First we define a continuous, 
strictly increasing, non- negative process, /3t, as 



/3t := / Zsds. 
Jo 

Next, we note that W^ may be dccompoascd as 



W,'=Py.W,' + Jl-plW,^, 



(C.l) 



where W{^ is a Brownian motion which is independent of W^ . Using equations 
(|5.18p and (|C.1| we derive 



N<Cr + ^ 



T^A 



ei^^Jz'sdW, 



+ e^^* 



ei^^Jz:dW' 



where Ci and C2 are constants. We will focus on bounding the first moment of 
the second stochastic integral. We have: 



Ie 

e 



^^* / ei^'y^dW, 



Ie 

e 


e-2A/e E 


( / e^^/S 


/ 


^.dVK,^) 


Pt 


Ie 


e-2A/^ E 


Ce^^'/'Zsds 


Pt 




Ie 

e 


e-2/5*/^ E 


' fe^P^/^dp, 


a] 




Ie 

e 


"e-2AAl (e^AA- iV 




^- 


"l-e-tA" 


<1. 

~ 2 













Then, by the Cauchy-Schwarz inequality, we see that 



=E 



^A 



e-'f^^JzldW. 



< 



V2- 



What remains is to bound the first moment of the other stochastic integral. 



A 



1 



=E 



ii/3* 



•}'^^JYsdw: 



Naively, one might try to use the Cauchy-Schwarz inequality in the following 
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manner 



A< ^w'E[e-2AA] Je 



s2/3=AZ,ds 



^ /e [e-2/5t/^] ^/e [- (e2^'/^) 



However, this approach does not work, since E [e^'^*/'^] — >■ oo as e ^> 0. Seeking 
a more refined approach of bounding A^ we note that 



1 -^A /".T/3=,/7rfW^ 



v^ 



— ^-^A(Z,-2)-^e-^^* / e^'5= 



crVe 



— Zs)ds 



,3/2 



Jo 



4A / .^/3. 



which can be derived by replacing t by s in equation (j2.4p . muhiplying by e^^'"^, 
integrating the resuh from to i and using Z^ = ZtZg — Zs{Zt — Zs) and 
/o G~'^^*~^''^^'^Zsds = e(l — e~^*/^). From the equation above, we see that 



A< 



a^e 



E 



'^'^'IZ, 



t-z\ 



ay/l 



E 



4a 



,kPs 



Zs)ds 



ae 



3/2 



E 



c'^^ZsiZt - Zs)ds 



(C.2) 



At this point, we need the moment generating function of {Zt,l3t). From 
we have 



E [e-^^*-^A] = e 



■■K,0(|lx,^i{t)-z^px,^L{t) 



(C.3) 



_2 
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(7+K)i/2 



Acr2 (eTt _ 1) + ^ _ K + e''*(7 + n) 



i^xAt) 



A (7 + « + eT*(7 - At)) + 2/x (e^' - 1) 
Acr2 (eT* - 1) + 7 - K + eT*(7 + k) ' 



7 = v K^ + 2cr2/i. 



Now, let us focus on the first term in equation (jC.2p . Using Cauchy-Schwarz, 
we have 



1 
uy/e 



E 



1 



CTVe 



e-^'/^|Z, ~z\ <^Je [e-2A/^] yipT^^. 
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From equation (jC.Sp one can verify 



E 



,-2AA 



— -K00o,2/j(t)-20O,2/e(t) ^ pCij ^ 



where C3 and C4 are constants. Since -^ep^l^ -^^ as e -> we see that 



=E 






<C5 



for some constant C5. 

We now turn out attention to the second term in equation (|C.2p . We have 



CTVe 



=E 



e-i(A-0.)(0_^^)rf, 



< — pE 
(TVe 



-T(A-'3»)Z,,d,S 



< 



< 






<C6 + 







■\(fit-Ps) 



_p-A/< 



d/3s 



k6I 

E 



E 



cr^/e 



t 



-iiPt-Mds 



-m-Mds 







K0 



=E 



Kb 



E 






,-7i0^-Mds 



Uo 



(^V^ Jo 
for some constant Cq. To bound the remaining integral we calculate 



ds, 



E 



■7(A-/9=) 



E 



E 



'7(/3t-A) 



-K9</'o,i/j(t-s)--Zs'0o,i/j(*-s) 



= exp (^ - K94>os/e{t - s) - k6'(/)^_o(s) - z?/'^_o(s)j (C.4) 
V'(s) := V'o,i/£(^-^)- 
Using the fact that (/'A,/i(^), ^a.^C*) > for any X,fJ.,t > 0, we see that 



E 



4 (A- A) 



< cxp (^- zi/j^Q^s)^ 



Hence 

nt 







=-|(A-/3=) 



ds < exp 

Jo 

=:/i+/2, 



[- zipij,ois))ds+ I exp 



(C.5) 
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where a € (1/2, 1). Using again ip\^^{t) > we deduce i/)j, q(s) > and therefore 



h < e". 
As for /i, we claim 

/i<C7exp(-C8e"-^), 
which is equivalent to showing there exists a constant C such that 

for all s e [0,t — e"]. To prove this claim, we note that for small e 

1 



(C.6) 



(C.7) 



(C.8) 



V'(s) =l/'0,l/e(^- S) 



aV2 /exP 



exp 



^(^--) 



^(t-.) 



v^ 



1 



where we have used 7 = -^/k^ + 2a'^ je ^ a^/2/ ^/e. A direct computation shows 
that "0(5) is a strictly decreasing in s with 



Now, we note that -^j; o(s) is given by 



i'^fii.s) 



2Kip{s) 



2k 



^2 (gKs _ ;l) ^(5) _,_ 2^6'^-' cr2 (gKs - 1) + 2Ke'"' / l}}{s) 



Since e**^ < e''*, and since, at worst, ■0(s) ~ cr^e""^, we conclude that there 
exists a constant C such that (jC.Sp . and therefore (|C.7p . hold. Hence, using 
equation (|C.5IIC.7)) . we have 



--A0t-tis) 



ds 



<C-je 



E 



-Cse° 



'i(A-/3,) 



+ e" 



ds 



E 



-i(/3t-/3.) 



ds 



This implies that there exists a constant Cg such that for any a E (1/2, 1) 



E 



4(0*-/3-) 



ds <Cn 



Having established a uniform bound on the first two terms in equation (|C.2[) . 
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we turn our attention toward the third and final term. For a G (1/2, 1) we have 



1 



f3/2 



E 



■^(A-Ps) 



Zs{Zt - Zs)ds 



< 



1 



(je3/2 
1 

1 



E 



E 



,-i(A-/3.) 



Zs\Zt - Zs\ds 



+ 



CTe3/2 

For the integral from to (t — e") we compute 
1 , 



/ ' e-^^^'~^^'^Zs\Zt^Zs\ds 
Jo 

f e-i(^^-^^^Zs\Zt-Zs\ds 



p3/2 



-¥fi^-p^)z,\Zt-ZMs 



< 



3/2 ' 







" / 






N 2] 


\ 


E 


sup Zs 

\0<s<t 


\Zt- 


~Zs\ 


)] 



ae 



- ^372^10'^ "' - '^12' 



MPt-Mds 



for some constants Cio, Cn and C12. For the integral from (i — e") to < we have 



1 



ae 



3/2 



E 



.-|(A-/3=) 



Zs\Zt- Zs\ds 



< 



1 



< 



ae3/2 

1 

cre3/2 

1 



sup l^t ~ Zs\ I e 

t-t^KsKt Jt-e 



ae 



1/2 



sup |Zt-Z,|e(l-e-^('-'°)/0 



t-£a<S<i 



sup |Zt - Zs 
t-e°'<s<t 



< Ci3e° 



for some constant C13. With this result, wc have established that for all a G 
(1/2, 1) there exists a constant, C, such that E|yt| < Ce"~^. 

D Numerical Computation of Option Prices 



The formulas (j4.3p and (|4.19p for Pnit, x, z) and Pi{t, x, z) cannot be evaluated 
analytically. Therefore, in order for these formulas to be useful, an efficient 
and reliable numerical integration scheme is needed. Unfortunately, numerical 
evaluation of the integral in (j4.3|) is notoriously difficult. And, the double and 
triple integrals that appear in (|4.19p are no easier to compute. In this section, we 
point out some of the difficulties associated with numerically evaluating these 
expressions, and show how these difficulties can be addressed. We begin by 
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establishing some notation. 

P'{t,x,z) ^ PH{t,x,z) + VePiit,x,z), 





27r Jr 


3-^^* (1 + x/i (K0/o(r,fc) +z7i(r,fc))) 






xGiT,k, 


z)/i(fc)dfc, 






g-rr 

= -^ (PoAt^ 2;, z) + K6'^/iPl,o(t, X, z) + zyfePi,i{t, x, 


^)), 


where 


we have defined 








Po,o(i,a;,z) : 


= [ e-''"^G{T,k,z)h{k)dk, 

JtS. 


(D.l) 




Pi,o(t,a;,2;) : 


= 1 e-''"iUT,k)G{T,k,z)h{k)dk, 


(D.2) 




Pi,i{t,x,z) : 


= f €-''"> h{T, k)G{T, k, Z)h{k)dk. 

Jr 


(D.3) 


As the 


v are written. (|D.l(). 


(|D.2|) and (|D.3p are general enoueh to accc 


imodate 



any European option. However, in order to make progress, we now specify an 
option payoff. We will limit ourself to the case of an European call, which has 
payoff /i(x) = (x—K)^. Extension to other European options is straightforward. 
We remind the reader that h{k) is the Fourier transform of the option payoff, 
expressed as a function of g = r{T — t) + log(x). For the case of the European 
call, we have: 

Jr i/s /s 

We note that (|D.4p will not converge unless the imaginary part of k is greater 
than one. Thus, we decompose k into its real and imaginary parts, and impose 
the following condition on the imaginary part of k. 

k, > 1. (D.5) 

When we integrate over k in ()D.1|) . (jD.2p and (jD.3|) . we hold ki > 1 fixed, and 
integrate fc^ over R. 

Numerical Evaluation of PqA^^x^z) 

We rewrite (jD.ip here, explicitly using expressions (j4.7p and (|D.4p for G{t, k, z) 
and h{k) respectively. 

PoAt.x,z) = / e-*'^«e^(^''=)+^^(^''^-)— -— dfc,. (D.6) 

./Iff ^"^' f^ 
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In order for any numerical integration scheme to work, we must verify the con- 
tinuity of the integrand in (|D.6|) . First, by (jD.Sp . the poles at fc = and k = i 
are avoided. The only other worrisome term in the integrand of (|D.6p is e'-^'^'^''^', 
which may be discontinuous due to the presence of the log in C{t, fc). 

We recall that any ( G <C can be represented in polar notation as C = 
rexp{i9), where 9 £ [— 7r,7r). In this notation, logC = logr + i9. Now, suppose 
we have a map C(^r-) : R — >■ C We see that whenever C(fcr) crosses the negative 
real axis, log^(fcr) will be discontinuous (due to 9 jumping from —n to tt or 
from TT to — tt). Thus, in order for logC(fcr) to be continuous, we must ensure 
that ({kr) does not cross the negative real axis. 

We now return our attention to C{T,k). We note that C{T,k) has two 
algebraically equivalent representations, (|4.8p and the following representation: 



C{T,k) = —{{K + pika-d{k))T-2logC{T,k)), (D.7) 



a^ 



It turns out that, under most reasonable conditions, C(r, fc) does not cross the 
negative real axis |17j . As such, as one integrates over fc^, no discontinuities will 
arise from the log C{t, k) which appears in (|D.7[) . Therefore, if we use expression 
(|D.7p when evaluating (jD.6p . the integrand will be continuous. 

Numerical Evaluation of Pi i(t,x, 2) and Pi^Q{t,x,z) 

The integrands in (|D.3|) and (jD.2p arc identical to that of (|D.1[) . except for 
the additional factor of /i(t, k). Using equation (|4.2ip for /i(t, k) we have the 
following expression for Pi^i(t, x, z): 



Pi^i{t,x,z) 

\Jo ) ik-k^ 






■'^'6(S, k)e^{--^^^)^C{rM) + zD(.,\.) JL^^dkrds. (D.9) 

ik — k^ 



Similarly: 



Plfi{t,X,z) 

f-T ft f fy-l+ik 

-^kqf^^^^^f.yA{t,k,s)+C{r,k)+zD{rM) dkrdsdt.{B. 10) 

ik — fc^ 

We already know, from our analysis of Po,o{t, x, z), how to deal with the log in 
C{t, fc). It turns out that the log in A(t, fc, s) can be dealt with in a similar man- 
ner. Consider the following representation for A(r, fc,s), which is algebraically 
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equivalent to expression (|4.22p : 

X (d(fc)(T - s) + logC(T, k) - logC(s, fc)) 
+d(fc)(r-s), (D.ll) 

where ({t, k) is defined in (|D.8I) . As expressed in (jD.lip . A{t, k, s) is, under most 
reasonable conditions, a continuous function of fc^. Thus, if we use (jD.lip when 
numerically evaluating (|D.9|) and (|D.10p . their integrands will be continuous. 

Transforming the Domain of Integration 

Aside from using equations (|D.7P and (jD.lip for C{T,k) and A{T,k,s), there 
are a few other tricks we can use to facilitate the numerical evaluation of (ID.6P , 
(|D.10[) . and (jD.9p . Denote by Io{k) and Ii{k,s) the integrands appearing in 
(pel) . (IE3) and (iDlOl) . 



Po.o = / Io{k)dkr, 

Pi I — Ii{k,s)dkrds, 

Jo Jr 

Pi^o = / / / Ii{k,s)dkrdsdt. 
Jo Jo Jm. 



First, we note that the real and imaginary parts of /o(^') and /i(fc, s) are even 
and odd functions of kr respectively. As such, instead of integrating in kr over 
R, we can integrate in kr over M+, drop the imaginary part, and multiply the 
result by 2. 

Second, numerically integrating in kr over M_|_ requires that one arbitrarily 
truncate the integral at some kcutoff- Rather than doing this, we can make the 
following variable transformation, suggested by J14j : 

_ -logM 

Kr — r< ^ 



Coo := — —{z + Ker). (D.12) 

(T 



Then, for some arbitrary I(k) we have 



logu 



I{k)dkr = I I J= + ih -T^du. 

Jo \ ^oo / uL-oo 

Thus, we avoid having to establish a cutoff value, kcutoff (and avoid the error 
that comes along with doing so). 

Finally, evaluating (jD.10|) requires that one integrates over the triangular 



36 



region parameterized by < s < t < t. Unfortunately, most numerical inte- 
gration packages only facilitate integration over a rectangular region. We can 
overcome this difficulty by performing the following transformation of variables: 

s ~ tv, 
ds — tdv. 

Then, for some arbitrary I{s) we have 

/ / I{s)dsdt = I{tv)tdvdt. 

Jo Jo Ja Jo 

Pulling everything together we obtain: 

Po.o = 2Re f h f — 5^ + ik^ -^du, 

Jo V '-'oo / UC-oo 

Pij = 2Re I f h 1 !?^" +iKs] -^duds, 



(D.13) 



Jo .^0 \ ^oo / ^^CXD 

Pi,o = 2Re h ( ^ ^^^ + iki,tv\ -^dudvdt, 

Jo Jo Jo \ ^oo / uGoo 

where Coo is given by (jD.121) . These three changes allow one to efficiently and 
accurately numerically evaluate (|D.6|) . (|D.9|) and (jD.lOp . 

Numerical tests show that for strikes ranging from 0.5 to 1.5 the spot price, 
and for expirations ranging from 3 months to 3 years, it takes roughly 100 times 
longer to calculate a volatility surface using the multi-scale model than it does 
to calculate the same surface using the Heston model. 
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